

**********************************************************
*** Welfare simulation: 20x20 grid of kappa and sigma  ***
**********************************************************

* Load data, based on 4950 grid

capture ssc install matsave  

cd "C:\Users\hy65byfe\Desktop\smerge_0712"

* cd "R:\WSV2\TBu_AKe\Spatial_NEW"

use spatial_2_collapse_s, replace    

cd "R:\WSV2\TBu_AKe\Spatial_NEW\Welfare"

log using welfare_simulation_20x20, replace 

keep if year == 2011 // 

sort cell 

gsort -bunching +cell 
gen bcell = _n if bunching == 1 // identify B

gsort -restricted +cell 
gen rcell = _n if restricted == 1 // identify R

*** policy variables ***

scalar E_bar = 1000/220 // 

scalar K = (51.7*0.59)/220

scalar kappa = E_bar - `=K' // at compliance line, in days. 

scalar sigma = (47*0.59)/220 // at compliance line, in days


*** parameters ****
	
scalar beta = 0 // by assumption 

scalar gamma = sigma // from delta = 0: scalar delta = 1/(-tan((90*(_pi/180)))) 

gen E_bin = (k_bin/220) + a_bin*`=sigma' // 
 
replace eei_bin = ((E_bin*220)/(51.7+47*a_bin))* 100 //  


*****************************************************
*****************************************************
*** 1) preference parameters


gen a_0 = 		a_bin 						// 

gen e_0 = 		E_bar 	- E_bin 			//  

gen alpha = 	a_0 	+ beta 	+ e_0*gamma // beta drops out 

gen epsilon = 	e_0 	+ a_0*gamma 		//  


*****************************************************
*****************************************************
*** 2) predicted adjustment 

gen e_gap = kappa - e_0 - sigma*a_0 // 

gen a_gap = (sigma - gamma)*e_gap    // drops to zero if sigma = gamma. 

gen e_star = 	e_gap + e_0 					//  

gen a_star = 	a_gap + a_0 					//   


*****************************************************
*****************************************************
*** 3) evaluate Delta V. 

gen V_gap = 0.5*(e_star - e_0)^2  // 

gen E_star = E_bar - e_star  


*****************************************************
*****************************************************
*** 4.a) sum up weights (i.e. sales volume)

	*** aggregate 
sort year restricted 

by year, sort: egen w_sum_T =sum(cell_year_sales)  //  

by year restricted, sort: egen w_sum_R =sum(cell_year_sales) if restricted== 1 // by year

gen R_ratio = w_sum_R/w_sum_T

by year, sort: egen p_sum_T =sum(cell_year_count) 

by year restricted, sort: egen p_sum_R =sum(cell_year_count) if restricted == 1 // by year
	
gen P_ratio = p_sum_R/p_sum_T
	
*****************************************************
*****************************************************
*** 4.b) sum up Delta V 

	*** at cell level 

gen V_sum_cell = V_gap*cell_year_sales 

	*** aggregate 
sort year restricted 

by year restricted, sort: egen V_sum_R =sum(V_sum_cell) if restricted== 1 // by year

 
sum V_gap if restricted == 1 & year == 2011 // 

sum V_gap V_sum_cell V_sum_R if restricted == 1 & year == 2011 

codebook V_sum_R // 

tab cell a_bin if V_gap < 0 & restricted == 1 & year == 2011 // 
 
*****************************************************
*****************************************************
*** 5.a) sum up changes in E  and a 

gen E_gap = E_bin -  E_star // 

sum E_gap a_gap if restricted == 1 & year == 2011 // 

gen E_gap_cell = E_gap*cell_year_sales 

sort year restricted 

by year restricted, sort: egen E_sum_R =sum(E_gap_cell) if restricted== 1 // by year

gen E_mean_R = E_sum_R/w_sum_R // mean distance among restricted consumers 
		
gen a_gap_cell = a_gap*cell_year_sales 
		
sort year restricted 

by year restricted, sort: egen a_sum_R =sum(a_gap_cell) if restricted== 1 // by year

gen a_mean_R = a_sum_R/w_sum_R // mean distance among restricted consumers 


*****************************************************
*****************************************************
*** 5.b) construct variance 

 	*** at cell level for E 

gen D_sum_cell = cell_year_sales*(E_gap)^2 // squared sum of deviations, times sales volume in cell. 

	*** aggregate in R 
sort year restricted 

by year restricted, sort: egen D_sum_R =sum(D_sum_cell) if restricted== 1 // sum over R. 

	*** mean of the above. 
gen D_mean_R = D_sum_R / w_sum_R // mean of squared deviations

 
sum D_sum_cell D_mean_R if restricted == 1 & year == 2011 // 
	
gen DA_sum_cell = cell_year_sales*(a_gap)^2 // squared sum of deviations, times sales volume in cell. 

	*** aggregate in R 
sort year restricted 

by year restricted, sort: egen DA_sum_R =sum(DA_sum_cell) if restricted== 1 // sum over R. 

	*** mean of the above. 
gen DA_mean_R = DA_sum_R / w_sum_R // mean of squared deviations

sum DA_sum_cell DA_mean_R if restricted == 1 & year == 2011 // 

*****************************************************
*****************************************************
*** 5.c) construct coefficient of variation 

 	*** at cell level 

gen cv_R = sqrt(D_mean_R)/E_mean_R // std. deviation / mean. 


sum cv_R D_mean_R E_mean_R if restricted == 1 & year == 2011 // 


*****************************************************
*****************************************************
*** 6) utility ratio  


gen U_ratio = V_sum_R / E_sum_R 

sum w_sum_R E_sum_R E_mean_R D_mean_R V_sum_R U_ratio a_sum_R a_mean_R DA_mean_R cv_R restricted if year == 2011 //  

* restore  

	*****************************************************
	** 				Simulation over sigma and kappa    **
	*****************************************************
	
*** Matrices for outcomes 
	mat define V = I(20) // utility ratio 
	mat define E = I(20) // total E savings   
	mat define W = I(20) // fraction of sales 
	mat define P = I(20) // fraction of products 

	
forvalues i = 1(1)20 { 	
	
	scalar kappa = 3+`i'*0.1-0.1 // vary kappa from 3 to 5 in 0.01 increments
	local k = kappa 
	
	forvalues j = 1(1)20 {
	
	scalar sigma = 0+`j'*0.01-0.01 // vary sigma from 0 to 0.20 in increments of 0.01 	
	local s = sigma 
	
		*** drop variables from previous loop 
	drop e_gap-U_ratio // reset to overwrite 
	
	gen restricted_`i'_`j' = 0 
	replace restricted_`i'_`j' = 1 if e_0 < kappa - sigma*a_0
	
	**** 
	*****************************************************
	*****************************************************
	*** 2) predicted adjustment 


	gen e_gap =  	( (kappa - e_0 - sigma*a_0)	/	(1 - 2*sigma*gamma + sigma^2) )   * ( 1 - sigma*gamma) // equivalent to E_gap. 
	
	replace e_gap = 0 if e_gap < 0 // remove values with negative adjustment 
	
	gen a_gap  =  -1* ( ( (gamma - sigma)* (kappa - e_0 - sigma*a_0)  )   /   (1 - 2*sigma*gamma + sigma^2) )

	gen  e_star = 	e_gap + e_0 					// 

	gen  a_star = 	a_gap + a_0 					//    
 
	*****************************************************
	*****************************************************
	*** 3) evaluate Delta V. 

	gen  V_gap = 0.5*(1-gamma^2)/(1-2*sigma*gamma+sigma^2)*(e_star - e_0)^2  //  
	
	gen  E_star = E_bar - e_star  


	*****************************************************
	*****************************************************
	*** 4.a) sum up weights (i.e. sales volume)

	*** added 05/02/21 to reflect fraction of sales  
	
	sort year restricted_`i'_`j' 
	
	by year, sort: egen w_sum_T =sum(cell_year_sales)  // 

	by year restricted_`i'_`j', sort: egen w_sum_R =sum(cell_year_sales) if restricted_`i'_`j'== 1 // by year

	gen R_ratio = w_sum_R/w_sum_T

	by year, sort: egen p_sum_T =sum(cell_year_count) 

	by year restricted_`i'_`j', sort: egen p_sum_R =sum(cell_year_count) if restricted_`i'_`j' == 1 // by year
	
	gen P_ratio = p_sum_R/p_sum_T


	*****************************************************
	*****************************************************
	*** 4.b) sum up Delta V 

	*** at cell level 

	gen V_sum_cell = V_gap*cell_year_sales 

	*** aggregate 
	
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen V_sum_R =sum(V_sum_cell) if restricted_`i'_`j' == 1 // by year

 
	sum V_gap if restricted_`i'_`j' == 1 & year == 2011 // 

	sum V_gap V_sum_cell V_sum_R if restricted_`i'_`j' == 1 & year == 2011 

	codebook V_sum_R // 

	tab cell a_bin if V_gap < 0 & restricted_`i'_`j' == 1 & year == 2011 // 
 


	*****************************************************
	*****************************************************
	*** 5.a) sum up changes in E  and a 

	* gen E_gap = E_bin -  E_star // 600 drops out: e_gap = 600 - e_star - 600 + e_0. 
	
	gen E_gap = e_star - e_0 

	sum E_gap a_gap if restricted_`i'_`j' == 1 & year == 2011 // 
	gen E_gap_cell = E_gap*cell_year_sales 

	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen E_sum_R =sum(E_gap_cell) if restricted_`i'_`j'== 1 // by year

	gen E_mean_R = E_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 
		
	gen a_gap_cell = a_gap*cell_year_sales 
		
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen a_sum_R =sum(a_gap_cell) if restricted_`i'_`j'== 1 // by year

	gen a_mean_R = a_sum_R/w_sum_R // mean distance among restricted_`i'_`j' consumers. 


	*****************************************************
	*****************************************************
	*** 5.b) construct variance 

 	*** at cell level for E 

	gen D_sum_cell = cell_year_sales*(E_gap)^2 // squared sum of deviations, times sales volume in cell. 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen D_sum_R =sum(D_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen D_mean_R = D_sum_R / w_sum_R // mean of squared deviations

 
	sum D_sum_cell D_mean_R if restricted_`i'_`j' == 1 & year == 2011 // 

	
	gen DA_sum_cell = cell_year_sales*(a_gap)^2 // squared sum of deviations, times sales volume in cell. 

	*** aggregate in R 
	sort year restricted_`i'_`j' 

	by year restricted_`i'_`j', sort: egen DA_sum_R =sum(DA_sum_cell) if restricted_`i'_`j'== 1 // sum over R. 

	*** mean of the above. 
	gen DA_mean_R = DA_sum_R / w_sum_R // mean of squared deviations

 
	sum DA_sum_cell DA_mean_R if restricted_`i'_`j' == 1 & year == 2011 //  


	*****************************************************
	*****************************************************
	*** 5.c) construct coefficient of variation 

 	*** at cell level 

	gen cv_R = sqrt(D_mean_R)/E_mean_R // std. deviation / mean. 


	sum cv_R D_mean_R E_mean_R if restricted_`i'_`j' == 1 & year == 2011 //  


	*****************************************************
	*****************************************************
	*** 6) utility ratio  


	gen U_ratio = V_sum_R / E_sum_R 

	sum w_sum_R E_sum_R E_mean_R D_mean_R V_sum_R U_ratio a_sum_R a_mean_R DA_mean_R cv_R restricted_`i'_`j' if year == 2011 //  

* restore 
	
	*****************************************************
	*****************************************************
	*** 7) transform into matrix  
	
	quietly sum E_sum_R if year == 2011 // outcome: savings 
	scalar e_`i'_`j' = r(mean) 
		
	quietly sum U_ratio if year == 2011 // outcome: v_loss per kwh 
	scalar v_`i'_`j' = r(mean)
	
	quietly sum R_ratio if year == 2011 // outcome: share of sales under restriction
	scalar w_`i'_`j' = r(mean)
	
	quietly sum P_ratio if year == 2011 // outcome: share of products under restriction
	scalar p_`i'_`j' = r(mean)
	
	** run through combinations and fill in outcomes **
			mat V[`i',`j']=(v_`i'_`j')
			mat E[`i',`j']=(e_`i'_`j')
			mat W[`i',`j']=(w_`i'_`j')
			mat P[`i',`j']=(p_`i'_`j')
	} 	
	
}	
	
	
*** matsave for each outcome ***

matsave V, replace saving
matsave E, replace saving
matsave W, replace saving
matsave P, replace saving

/* 
svmat V, names (v)
svmat E, names (e)
svmat W, names (w)
svmat P, names (p)

save welfare_3D_database 

*/ 

log close 
exit 
clear 
